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ABSTRACT 

We present a study of numerical effects in dissipationless cosmological simulations. The 
numerical effects are evaluated and studied by comparing results of a series of 64 3 - 
I I particle simulations of varying force resolution and number of time steps, performed 

. using three of the TV-body techniques currently used for cosmological simulations: the 

' Particle Mesh (PM), the Adaptive Particle-Particle Particle- Mesh (AP 3 M), and the 

I/" - ) I newer Adaptive Refinement Tree (ART) code. This study can therefore be interesting 

£NJ ■ both as an analysis of numerical effects and as a systematic comparison of different 

CN ' codes. 

We find that the AP 3 M and the ART codes produce similar results, given that 

ON . 

convergence is reached within the code type. We also find that numerical effects may 
affect the high-resolution simulations in ways that have not been discussed before. In 
r"| . particular, our study revealed the presence of two-body scattering, effects of which can 

(*S be greatly amplified by inaccuracies of time integration. This process appears to affect 

' \ the correlation function of matter, mass function and inner density of dark matter halos 

and other statistics at scales much larger than the force resolution, although different 
' statistics may be affected in different fashion. We discuss the conditions at which 

\ strong two-body scattering is possible and discuss the choice of the force resolution 

■ and integration time step. Furthermore, we discuss recent claims that simulations with 

, £^ ' force softening smaller than the mean interparticle separation are not trustworthy and 

, argue that this claim is incorrect in general and applies only to the phase-sensitive 

5_i ■ statistics. Our conclusion is that, depending on the choice of mass and force resolution 

and integration time step, a force resolution as small as 0.01 of the mean interparticle 
separation can be justified. 
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1 INTRODUCTION 

Dissipationless cosmological iV-body simulations are cur- 
rently the tool of choice for following the evolution of 
cold dark matter (CDM) into the highly nonlinear regime. 
For the widest range of plausible dark matter (DM) can- 
didates (from axions of mass ~ 10" 6 - 10" 3 eV to ~ 
10 21 - 10 25 eV wimpzillas; see, e.g., Kolb, Chung & Ri- 
otto 1998; Roszkowski 1999), their expected number den- 
sity is n D M ~ (10 52 - 10 83 )Q h 2 Mpc -3 . iV-body simula- 
tions numerically solve the A^-body problem: given initial 
positions and velocities for N pointlike massive objects, 
the simulations predict the particle positions and veloci- 
ties at any subsequent time. Current Af-body simulations 



are capable of following the evolution of < 10 particles, 
far short of the expected number of DM particles. There- 
fore, the correct approach to modelling dark matter evolu- 
tion in a cosmologically representative volume, is to use the 
Vlasov equation (collisionless Boltzmann equation) coupled 
with the Poisson equation and complemented by appropriate 
boundary conditions. However, a full-scale modelling of 6D 
distribution functions with reasonable spatial resolution is 
extremely challenging computationally. The alternative ap- 
proach adopted in most cosmological simulations is to split 
the initial phase space into small volume elements and fol- 
low evolution of these elements using iV-body techniques. 
Each volume element can thus be thought of as an A^-body 
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particle, which moves with a flow and which has some shape 
(for example, a box or a sphere) or is simply point-like. Two 
particles are assumed to interact gravitationally only if they 
are separated by a distance > e, where e is the smoothing 
scale, often referred to as the force resolution. It is clear 
that the "size" of a particle is determined by the eulerian 
spatial size of the initial phase space volume element or by 
the smoothing scale, whichever is smaller. 

While this approach seems logical and reasonable and is 
expected to provide approximate solution to a complicated 
problem, questions of its limitations may be raised. For ex- 
ample, while particle shape is usually considered to be rigid 
(fixed by a specific form of the Green function or the shape 
of the interparticle force), in eulerian space the shape of the 
initially cubic phase space volume element can be expected 
to be stretched as the element moves towards higher den- 
sity regions. Its volume can also change so as to preserve 
the phase space density. Furthermore, under certain condi- 
tions the iV-body systems may exhibit scattering, which is 
undesirable when one models a purely collisionless system. 
This may occur in cosmological simulations if the "size" of 
particles is much smaller than the eulerian spatial size of 
the phase space element they are supposed to represent. 

These effects may influence the accuracy of the simula- 
tions and lead to spurious results. Nevertheless, surprisingly 
little attention has been given to studies of such limitations 
in the cosmological simulations. As the resolution of simula- 
tions improves and the range of their applications broaden, 
it becomes increasingly important to address these issues. 
Indeed, during the past decade the force resolution of the 
simulations has improved by a factor of ~ 100 — 1000, while 
(with rare exceptions) the mass resolution has improved only 
by a factor of ~ 10. Also, modern high-resolution codes fol- 
low evolution of cosmological systems for many dynamical 
time-scales. In this regime the accuracy of the force esti- 
mates may be less important than the stability of the overall 
solution. These issues are usually not addressed when tests 
of codes are presented. 

Recently, in a series of papers Melott and collaborators 
(Kuhlman, Melott & Shandarin 1996; Melott et al. 1997; 
Splinter et al. 1998) raised the issue of the balance between 
the force and the mass resolution. While we disagree with 
their main conclusion that the force resolution should be 
larger than the mean interparticle separation (see § 3 and 
§ 5), we agree that the issue is important. 

There is also a common misconception related to the 
adaptive mesh refinement approach in cosmological simula- 
tions and other algorithms that integrate equations of parti- 
cle motion in comoving coordinates. The common criticism 
(e.g., Splinter et al. 1998 and references therein) is that these 
algorithms attempt to resolve scales unresolved in the ini- 
tial conditions (the scales below approximately half of the 
Nyquist frequency). However, the goal of increasing the res- 
olution in comoving coordinates is not to resolve the waves 
not present in the initial conditions but rather to properly 
follow all of the waves initially present. 

When followed in comoving coordinates, gravitational 
instability leads to the separation of structures from the 
Hubble flow and collapse, resulting in the transfer of power 
to higher wavenumbers. If the force resolution is fixed in 
comoving coordinates at the Nyquist frequency of the ini- 
tial conditions, this transfer cannot be modelled properly 



for all waves. Moreover, the size of structures that have col- 
lapsed and virialized stays fixed in the proper coordinates 
and decreases in comoving coordinates. When the comoving 
size of such objects becomes smaller than the resolution of a 
fixed grid simulation, their subsequent evolution and inter- 
nal properties will be modelled incorrectly. Adaptive mesh 
refinement algorithms address this problem by increasing 
the force resolution locally to follow the evolution of col- 
lapsing and virialized density peaks as their size becomes 
less than the resolution of the original grid. Other algorithms 
can achieve the same result by varying the comoving force 
softening with time. 

Our motivation for the present study is twofold. Our 
first goal is to elaborate on the issue of spurious numer- 
ical effects. Namely, we study the effects of the balance of 
force and mass resolutions and of time integration details on 
statistics commonly used in analyses of cosmological simu- 
lations. The balance of force and mass resolutions should 
be studied by varying both of the resolutions. In this study, 
however, we will keep the mass resolution fixed and vary the 
force resolution instead. While this may not reveal all of the 
artificial effects^, this allows to study the spurious two-body 
scattering present when the smoothing scale is set to be too 
small. It also allows us to study how the limited force resolu- 
tion affects various statistical properties of the dark matter 
distribution. Our second goal is to compare results produced 
using two of the currently employed high-resolution A-body 
codes: AP 3 M (Couchman 1991) and ART (Kravtsov, Klypin 
& Khokhlov 1997). The comparison is of interest due to a rel- 
ative novelty of the latter technique and some disagreement 
between the results concerning the details of the central den- 
sity distribution in DM halos obtained using different codes 
(e.g., Moore 1994; Navarro, Frenk & White 1997; Kravtsov 
et al. 1998; Moore et al. 1999). This study can also be inter- 
esting as an independent test of the widely used AP 3 M code, 
for which virtually no systematic tests have been published 
to date. 

The paper is organized as follows. In §2 we describe the 
A-body algorithms used in our study and describe the nu- 
merical simulations performed . In §3 we compare the dark 
matter distribution simulated using different codes with dif- 
ferent force resolutions and time steps. In §4 we use these 
simulations to compare the properties and distribution of 
dark matter halos (dense virialized systems). In §5 we sum- 
marize the results of the code comparisons, discuss the ef- 
fects of resolution and time step on the commonly used 
statistics, and present our conclusions. 



2 COSMOLOGICAL A-BODY SIMULATIONS 

In this study we will use and compare three different N- 
body algorithms: Particle-Mesh (PM) algorithm (Hockney 
& Eastwood 1981), adaptive Particle-Particle Particle- Mesh 
algorithm (AP 3 M; Couchman 1991), and Adaptive Refine- 
ment Tree algorithm (ART; Kravtsov et al. 1997). The PM 
algorithm was first used for cosmological simulations by 
Doroshkevich et al. (1980), Efstathiou & Eastwood (1981), 

* For example, fixed mass resolution does not allow us to study 
the effects of the rigid particle shape mentioned above. 
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and Klypin & Shandarin (1983). The algorithm makes use 
of the fast Fourier transforms to solve the Poisson equation 
on a uniform grid and uses interpolation and numerical dif- 
ferentiation to obtain the force that acts on each particle. 

The solution is limited by the number of particles (mass 
resolution) and by the size of the grid cell which defines force 
resolution. The exact shape of the resulting force depends 
on the specific form of the Green function and interpola- 
tion used to get the force. The technqiue is attractive due 
to its simplicity and the fact that it is numerically very ro- 
bust. Highly efficient implementations have been developed 
and used during the past decade. The technique is described 
in detail by Hockney & Eastwood (1981) and we refer the 
reader to this book for further details. The specific imple- 
mentations used in our study are those of the AP 3 M code 
and the ART code. The PM simulations presented here have 
been run using the publicly available AP 3 M code with the 
particle-particle part switched off and by the ART code with 
the mesh refinement block switched off. In the remainder of 
this section we will describe the AP 3 M and ART algorithms 
and the specifics of our test simulations. 

2.1 AP 3 M Code 

Particle-Particle-Particle-Mesh (P 3 M) codes (Hockney et al. 
1973; Hockney & Eastwood 1981) express the inter-particle 
force as a sum of a short range force (computed by direct 
particle-particle pair force summation) and the smoothly 
varying part (approximated by the particle-mesh force cal- 
culation). One of the major problems for these codes is the 
correct splitting of the force into a short-range and a long- 
range part. The grid method (PM) is only able to produce 
reliable inter particle forces down to a minimum of at least 
two grid cells. For smaller separations the force can no longer 
be represented on the grid and therefore one must introduce 
a cut-off radius r e (larger than two grid cells !) where for 
r < r e the force should smoothly go to zero. The parameter 
r e defines the chaining-mesh and for distances smaller than 
this cutoff radius r e a contribution from the direct particle- 
particle (PP) summation needs to be added to the total force 
acting on each particle. Again this PP force should smoothly 
go to zero for very small distances in order to avoid unphys- 
ical particle-particle scattering. This cutoff of the PP force 
determines the overall force resolution of a P 3 M code. 

The most widely used version of this algorithm is cur- 
rently the adaptive P 3 M (AP 3 M) code of Couchman (1991). 
The smoothing of the force in this code is connected to a S2 
sphere, as described in Hockney & Eastwood (1981). The 
particles are treated as density spheres with a profile 

f fl-rV for r<e/2 

1 0, otherwise 

where e is the softening parameter. For distances greater 
than e the particles are treated as point masses interacting 
according to the newtonian 1/r 2 law, whereas for smaller 
separations the effective shape of the S2 sphere influences 
and modifies the force law such that the interaction drops 
down to zero as r — > 0. 

When splitting the force into short- and long-range com- 
ponents, one has to use two softening parameters: one which 



is directly connected to the cut-off radius r e for the PM force 
and therefore tells us where to match the PM and PP part, 
and another which determines the overall force resolution 
(softening scale for PP force). The PP force is truncated 
at both the very low separations and at r > r e where the 
force can be calculated using the mesh based PM method. 
The AP 3 M code uses a cut-off radius r e for the long-range 
force of approximately 2.4 PM mesh cells, and this leads to 
the softening parameter of epm = 1.3r e » 3.1 (cf. Hockney 
& Eastwood 1981). The softening e of the PP force deter- 
mines the overall force resolution; for the AP 3 M simulations 
presented in this paper the softening scales are given in Ta- 
ble |l|. When setting the overall softening parameter to a 
value greater than 3.5 the code runs as a pure PM mesh 
based code, because the softening of the PP force is greater 
than that of the PM part. 

Unfortunately, the particle-particle summation which 
allows one to achive sub-grid resolutions and thereby makes 
the P 3 M algorithm attractive, is also the method's main 
drawback. The PP calculation must search for neighbors 
out to roughly two mesh spacings to properly augment the 
PM force. This becomes increasingly expensive as cluster- 
ing develops and particles start to clump together (within 
the chaining-mesh cells). The adaptive P 3 M algorithm reme- 
dies this by covering the high-density, most computationally 
expensive regions with refinement grids. Within the refine- 
ments the direct sum is replaced by a further local P 3 M 
calculation with isolated boundary conditions performed on 
a finer refinement grid (PM mesh and chaining-mesh are 
refined). The number of particles per refinement grid cell is 
smaller and so is the PP associated computations. The num- 
ber of grid cells per refinement depends on the total number 
of particles within that region, but is always a power of two. 
The criterion for placing a refinement depends only on the 
total number of particles inside a chaining-mesh cell. If this 
value exceeds a preselected threshold in a given region, the 
region is refined; for our runs we used a value of 50 particles. 
It is convenient to isolate patches which cover an exact cubic 
block of chaining-mesh cells. Recursively placed refinements 
are allowed, and in the simulations presented in this paper 
a maximum level of 3 was reached. 

The AP 3 M code integrates equations of particle motion 
using p = a 3 / 2 as a time variable (here a is the expansion 
factor; Efstathiou et al. 1985). A time-centered leapfrog in- 
tegration with constant time step Ap is used. This scheme 
leads to "large" steps in a at the beginning of the simula- 
tion which then get progressively smaller during the course 
of simulation. 

2.2 ART Code 

The Adaptive Refinement Tree code (ART; Kravtsov et 
al. 1997) reaches high force resolution by refining all high- 
density regions with an automated refinement algorithm. 
The refinements are recursive: the refined regions can also be 
refined, each subsequent refinement having half of the previ- 
ous level's cell size. This creates an hierarchy of refinement 
meshes of different resolutions covering regions of interest. 

The refinement is done cell-by-cell (individual cells can 
be refined or de-refined) and meshes are not constrained to 
have a rectangular (or any other) shape. This allows the code 
to refine the required regions in an efficient manner. The 
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criterion for refinement is the local overdensity of particles: 
in the simulations presented in this paper the code refined 
an individual cell only if the density of particles (smoothed 
with the cloud-in-cell scheme; Hockney & Eastwood 1981) 
was higher than nth = 5 particles. Therefore, all regions 
with overdensity higher than 6 = nth 2 3i /n, where n is the 
average number density of particles in the cube, were refined 
to the refinement level L. For the two ART simulations pre- 
sented here, n — 1/8. The Poisson equation on the hierarchy 
of meshes is solved first on the base grid and then on the 
subsequent refinement levels. On each refinement level the 
code obtains the potential by solving the Dirichlet boundary 
problem with boundary conditions provided by the already 
existing solution at the previous level. There is no particle- 
particle summation in the ART code and the actual force 
resolution is equal to » 2 cells of the finest refinement mesh 
covering a particular region. A detailed description of the 
code, its tests, and discussion of the force shape is given 
in Kravtsov et al. (1997). Note, however, that the present 
version of the code uses multiple time steps on different re- 
finement levels, as opposed to the constant time stepping in 
the original version of the code. The multiple time stepping 
scheme is described in some detail in Kravtsov et al. (1998; 
also see below). 

The refinement of the time integration mimics spatial 
refinement and the time step for each subsequent refinement 
level is two times smaller than the step on the previous level. 
Note, however, that particles on the same refinement level 
move with the same step. When a particle moves from one 
level to another, the time step changes and its position and 
velocity are interpolated to appropriate time moments. This 
interpolation is first-order accurate in time, whereas the rest 
of the integration is done with the second-order accurate 
time centered leap-frog scheme. All equations are integrated 
with the expansion factor a as a time variable and the global 
time step hierarchy is thus set by the step Aao at the zeroth 
level (uniform base grid). The step on level L is then Aol = 
Aa /2 L . 

The choice of an appropriate time step for a simulation 
is dictated by the peak force resolution. The number of time 
steps in our simulations is such that the rms displacement of 
particles during a single time-step is always less than 1/4 of 
a cell. No particles moves further than ~ 0.5 cells in a single 
time step, where the cell size and time step for particles 
located on the refinement level L are Axo/2 L and Aao/2 L , 
respectively. The value of Aao = 0.0015, used in run ART1 
(see Table 111) was determined in a convergence study using 
a set of 64 3 particle simulations described in Kravtsov et al. 
(1998). To study the effects of time step, we have also run a 
simulation with Aao = 0.003. 

The ART code integrates the equations of motion in 
comoving coordinates. Therefore, if a fixed grid is used to 
calculate the forces, the force resolution of the simulation 
degrades as oc a — (1 + z) _1 (see § 1). In order to pre- 
vent this and to preserve the initial resolution in physical 
coordinates in the simulations presented in this paper, the 
dynamic range between the start (zi = 87) and the end 
(z — 0) of the simulation should increase by (1 + Zi): i.e., for 
our simulations reach 128 x (1 + zC) — 11, 136. 

In the simulations presented in this paper, the peak res- 
olution is reached by creating a refinement hierarchy of five 
levels of refinement in addition to the base 128 3 uniform 



Table 1. Parameters of the numerical simulations. The corre- 
sponding force resolution in h~ x kpc is included in parenthesis. 
Note that PM1 and PM2 were simulated using the AP 3 M code, 
while PM3 and PM4 have been simulated using the ART code 
with mesh refinement switched off. The number of time steps for 
the ART simulations is given for the zeroth level; the effective 
number of time steps for particles on level L is x2 L , giving 21120 
steps for the maximum refinement level of the ART1 run. 



simulation 


softening 


dyn. range 


steps 


AP 3 M5 


0.06 


(7.0) 


2133 


8000 


AP 3 M1 


0.03 


(3.5) 


4267 


8000 


AP 3 M4 


0.03 


(3.5) 


4267 


2000 


AP 3 M2 


0.02 


(2.3) 


6400 


6000 


AP 3 M3 


0.015 


(1.8) 


8544 


6000 


ART1 


0.03125 


(3.7) 


4096 


660 


ART2 


0.03125 


(3.7) 


4096 


330 


PM1 


3.5 


(409) 


128 


6000 


PM2 


3.5 


(409) 


128 


1500 


PM3 


2 


(234) 


128 


1500 


PM4 


2 


(234) 


128 


750 



grid. However, the small number of particles in these simu- 
lations does not allow the code to reach the required target 
dynamic range of 11, 136, estimated above. 

2.3 Simulations 

Although we are not directly interested here in cosmological 
applications we decided to use a definite cosmological model: 
the cluster normalized, as = 0.67, standard cold dark mat- 
ter (SCDM) model. All simulations were run on a 128 3 mesh 
with 64 3 particles and started with the same random real- 
ization at z = 87. The adopted box size is 15/i _1 Mpc, which 
gives a mass resolution of 3.55 x 10 9 /i~ 1 M Q . The AP 3 M sim- 
ulations were carried out varying both force resolution and 
time steps. The two runs of the ART code differ only in the 
number of integration steps on the lowest-resolution of the 
uniform grid. 

For comparison we also ran the PM simulations. Both 
the AP 3 M and ART have an internal PM block and we ran 
two pairs of PM simulations using these different PM imple- 
mentations. This is done to compare implementations and 
to explore the effects of the time integration scheme on the 
results. In the case of AP 3 M (PM1 and PM2) both the adap- 
tive and the PP part of the P 3 M code have been switched 
off, while in the ART PM runs (PM3 and PM4) we have 
simply switched off mesh refinement. The parameters of the 
simulations are summarized in Table |l[ The force soften- 
ing is given in grid units and in /i _1 kpc (in brackets), and 
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Figure 1. Slice of 3 ft," 1 Mpc (top row) through the ART (left), AP 3 M (middle), and PM(right) simulation. In the bottom 
into the region marked in the top row. 



row wc zoom 



the number of steps of the ART simulations is presented for 
the lowest-resolution (the effective number on the highest- 
resolution is 32 times larger). 

Using this set of runs we can compare simulations with 
the same force resolution but different integration steps 
(e.g., AP 3 M1 with AP 3 M4, ART1 and ART2) amongst each 
other and simulations with the same number of integration 
steps but with varying force resolution (e.g., AP 3 M1 with 
AP 3 M5). Additionally we compare the different A^-body 
codes in order to quantify the deviations due to different 
(grid-based) methods to solve the equations of motion. 

It is important to keep in mind that the shape of the 
small-scale force is somewhat different in the codes used. 
Therefore, equal dynamic range does not correspond to the 
same physical resolution. The peak resolution of the ART 
code is ~ 2 cells of the highest level refinement, and so the 
actual force resolution is twice worse than the "formal" reso- 
lution given by the dynamic range. To make cross-code com- 
parison, we have performed the simulation AP 3 M5, which 
has approximately half the dynamic range of the ART runs 
and similar force resolution. 



3 PROPERTIES OF THE PARTICLE 
DISTRIBUTION 

3.1 Visual comparison 

A first inspection of the global distribution of particles in 
a 3/i _1 Mpc thick slice (Fig. [l], top row) shows that the dis- 
tributions are very similar. Even the much lower resolution 
PM1 simulation shows virtually the same global particle dis- 
tribution. In the bottom row of the figure we zoom into the 
region marked by rectangles in the upper panel. Here one 
can clearly see that the two high resolution runs produce 
many small-size dense halos, which are, however, slightly 
shifted in their positions. We attribute this shift to the cu- 
mulative phase errors due to the differences in the time in- 
tegration schemes of the type observed and discussed in a 
recent code comparison by Frenk et al. (1999; Santa-Barbara 
cluster comparison project). Most of these small clumps are 
absent in the PM run due to its poor resolution. The small 
scale density peaks do not collapse on scales smaller than the 
2 PM grid cells due to the sub-newtonian of self-gravitation 
at these scales (see § 3.3). From Figure [l] it can be seen that 
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Figure 2. Number of connections of the minimum spanning tree 
per bin of the connection lengths l con as a function of the con- 
nection length. The connection length is in units of the mean 
interparticle distance I, the bin size is dl con = 0.005 * I. 
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Figure 3. Left: Projections of 5% of all particles of the box along 
the z-axis. Right: Projections of all doublets and triplets found 
with 11 = 0.015. Note that the systems linked using this linking 
length have overdensity of 5 ^ 420000. 



only halos of size larger than approximately one grid cell 
collapse in the PM simulation (see, e.g., Klypin 1996). 

Our comparison can be contrasted with comparison by 
Suisalu & Saar (1996) (Fig.l and Fig. 2 in their paper). If 
compared on a large scale, the particle distributions in dif- 
ferent runs compare much better in our case than the sim- 



ulations compared in Suisalu & Saar (1996). This indicates 
that the differences they observed are due to the differences 
in their multigrid algorithm rather than to the two-body 
scattering, as argued in their paper. 

3.2 The minimal spanning tree 

To quantify the differences between the simulations we have 
calculated the minimal spanning tree (MST) of the particle 
distribution. The minimal spanning tree of any point distri- 
bution is a unique, well defined quantity which describes 
the clustering properties of the point process completely 
(e.g., Bhavsar & Splinter 1996, and references therein). The 
minimal spanning tree of n points contains N — 1 connec- 
tions. For the ART1, AP 3 M1 AP 3 M4, AP 3 M5 and a PM1 
simulations we show in Fig. H the number of connections 
N co „ of the tree per bin of the connection length (equal to 
0.005 * I) as a function of the length of the connection, l con - 
Here, T= (V^/iV) 1 / 3 denotes the mean inter-particle sepa- 
ration. Since the length of the connection is proportional to 
p _1//3 the probability distribution of connections (N co „/N) 
is equivalent to the density probability distribution in the 
simulation. 

In Fig. |i| the connection length distributions for ART 
and AP 3 M simulations are peaked at the same relative con- 
nection length (« 0.015 — 0.02/), whereas the PM simula- 
tion is peaked at a higher value (« 0.03 — 0.04Z). This in- 
dicates the ability of the ART and the AP 3 M codes to re- 
solve shorter scales and, therefore, to reach higher overden- 
sities. The position of the maximum depends only slightly 
on the resolution. In fact, increasing the resolution from PM 
to ART by a factor of 32 shifts the maximum by about a fac- 
tor of 2. This is probably because the differences affect only 
a very small fraction of the volume and dark matter parti- 
cles. Correspondingly, differences in resolutions of AP 3 M1 
and AP 3 M5 are too subtle to have a visible effect on the 
distribution. 

A time integration has a much more noticeable effect. 
An increase of the integration step in the run AP 3 M4 com- 
pared to AP 3 M1 and AP 3 M5 leads to a shift of the max- 
imum (and distribution) to higher values of l con - This is 
caused by inaccuracies in the integration of particle trajec- 
tories in high-density regions, as is evidenced, for example, 
by the halo density profiles in the AP 3 M4 run (see also re- 
sults below). 

The maxima of the AP 3 M simulations are slightly 
higher in amplitude than the maximum of the ART sim- 
ulation. This reflects the fact that the AP3M resolves forces 
uniformly (i.e., equally well in both low- and high-density re- 
gions) . The ART code by design reaches high resolution only 
in the high-density regions. Therefore, there are groups of a 
few particles located in low-density environments and thus 
not resolved by the ART code, which are however resolved 
by the AP 3 M. For example, we have found 3618 doublets 
and 558 triplets linked by = 0.015 (the maximum of the 
distribution) in the AP 3 M5 run, whereas only 2753 doublets 
and 466 triplets are found in the ART1 simulation at this 
linking lengths. 

Figure [] demonstrates that many of the missing dou- 
blets and triplets are indeed located in the low-density re- 
gions. The right column of this figure shows the projection 
of the distribution of all doublets and triplets found in the 
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Table 2. Density cross-correlation coefficient K = . 

C1IT2 



grid 

simulation 32 64 128 256 400 



ART1 <-> ART2 0.98 0.98 0.94 0.84 0.72 

AP 3 M1 <-» AP 3 M4 0.98 0.98 0.98 0.93 0.87 

PM1 «-> PM2 1.00 1.00 1.00 1.00 1.00 

PM3 «-» PM4 1.00 1.00 1.00 1.00 1.00 

AP 3 M2 «-» AP 3 M3 0.98 0.98 0.98 0.93 0.87 

AP 3 M1 <-> AP 3 M5 0.98 0.99 0.99 0.96 0.91 

PM2 <-> PM3 1.00 0.98 0.93 0.84 0.73 

ART1 <-> AP 3 M1 0.95 0.88 0.71 0.44 0.30 

ART1 <-» AP 3 M5 0.95 0.88 0.71 0.44 0.30 

ART1 <-> PM1 0.96 0.92 0.80 0.60 0.47 

AP 3 M1 <-> PM1 0.95 0.90 0.79 0.61 0.48 



AP 3 M5 (top) and ART1 (bottom) runs by the friends-of- 
friends algorithm with the linking length of 11 = 0.015. For 
comparison, we show in the left column of Fig. ^ a projec- 
tions of randomly selected 5% of all particles. While some 
doublets in the low-density regions are found in the ART 
run, all of them are gravitationally unbound (are chance su- 
perpositions). We find that in the AP 3 M runs most of such 
doublets and triplets are bound and are a part of a binary, 
triple, or higher multiplicity clusters. 

According to the sampling theorem, one needs at least 
~ 20 — 30 particles to resolve three-dimensional waves in 
the initial power spectrum. The presence of gravitationally 
bound clusters consisting of just a few particles is therefore 
artificial. 

3.3 Density Cross— Correlation Coefficient 

The density cross-correlation coefficient, 

K=^L, (2) 

0-1(72 

was introduced by Coles et al. (1993) in order to quantify 
similarities and differences between simulations of different 
cosmological models. Recently, Splinter et al. (1998) have 
adopted this statistic to compare simulations. Here, we fol- 
low the same approach and use this measure to quantify 
differences between simulations which have been carried out 
by different numerical algorithms or by the same algorithm 
but with different parameters of the simulation. 

To compute K, we have calculated the densities on a 
regular mesh using the triangular-shaped cloud (TSC; Hock- 
ney & Eastwood 1981) density assignment scheme and then 



used the resulting density field to compute (<5i<52) and vari- 
ances. We have varied the size of the grid in order to show 
the dependence of the cross-correlation on the smoothing 
scale of the density field. 

We summarize our results in Table ^: the first four 
rows present the cross-correlation coefficients between the 
runs with the same time integration scheme but different 
time steps. In the following two rows we present the cross- 
correlation coefficients for the AP 3 M runs with different 
force resolutions but the same integration step. The rest of 
the rows in the table present the cross-correlation coefficients 
for the runs with different time integration schemes as well 
as the cross-correlation coefficients between the AP 3 M1 and 
PM1 runs which have been simulated with the same time 
itegration scheme but with vastly different force resolutions. 

In all cases, it is obvious that the cross-correlation wors- 
ens for smaller smoothing scales (the larger density grid 
size). With smaller smoothing, smaller structures in the den- 
sity field are resolved. The degraded cross-correlation there- 
fore indicates that there are differences in locations and/or 
densities of these small-scale structures. It is also clear from 
the definition of K that this measure is particularly sensi- 
tive to the differences in the highest density regions. If we 
restrict the correlation analysis to a coarse grid, we smooth 
the particle distribution with a fairly large smoothing length 
and smear out the details and differences of the small-scale 
structure. 

The differences revealed by the cross-correlation coeffi- 
cient can arise either because the internal density distribu- 
tion of the structures is different in different runs or because 
the spatial locations of these structures are somewhat dif- 
ferent. Our analysis of halo profiles (See § 4.5) shows that 
differences in density in the same halos are small (except, 
of course for the PM runs) and have no significant effect 
on the cross-correlation. The degrading cross-correlation in 
high-resolution runs is thus due to the differences in the lo- 
cations of small-scale structure rather than to the differences 
in density. Indeed, one can readily see in the bottom row of 
Fig. [jjthat positions of small clumps in the ART and AP 3 M 
simulations often differ by ~ 100 — 300/i _1 kpc, the scale 
at which a significant decrease in K is observed. With such 
shifts in halo locations, the same halo may occupy different 
cells in the density grid which systematically reduces {8182) 
(calculated cell-by-cell) and results in a lower K. 

What causes the differences in halo positions? The rows 
1, 2, 5, and 6 of Table ^ indicate that both force resolution 
and time step have the same effect on K, both causing some 
phase errors. However, rows 7-9 show that the time integra- 
tion scheme causes much bigger phase differences than either 
force resolution or time step. This was indeed observed in 
the recent "Santa Barbara Cluster" code comparison project 
(Frenk et al. 1999). Different time integration schemes lead 
to a different accumulation of the phase errors. The mani- 
festation of these differences is certain "asynchronicity" be- 
tween simulations: the same phase error is accumulated at 
slightly different time moments. This results in shifts of halo 
positions when simulations are compared at the same time 
moment. 

We have indeed observed such asynchronicity in our 
simulations. Thus, for example, cross-correlation coefficient 
K between ART1 and AP 3 M5 on a 128 3 grid reaches a max- 
imum of 0.84 at a — 1.04 when AP 3 M5 is evolved further in 
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time (this can be compared to 0.71 at a = 1.0 in Table |2|) 
while keeping ART1 fixed at a=1.0. Similar effect is observed 
at z = 2: the maximum K is achieved when AP 3 M5 is ad- 
vanced forward in time by a factor of 1.02 in expansion fac- 
tor. Partly, this asynchronicity may be caused by the initial 
phase error introduced at the start as the AP 3 M simulations 
particles were advanced half a step (or by a factor of 1.03 in 
expansion factor) forward. The additional error accumulates 
during the time integration. 

While most halo properties and properties of the mat- 
ter distribution appear to be similar in the ART and the 
AP 3 M runs (see results below) , the differences between posi- 
tions of small-scale structures are much larger between these 
runs than any differences between runs simulated using the 
same code. This is clearly seen in the case of PM runs. All 
4 runs cross-correlate perfectly within the code type (PM1 
and PM2 were run using AP 3 M, while PM3 and PM4 were 
run using the ART code), but cross-correlate rather poorly 
when different code simulations are compared. In the lat- 
ter case we observed a decrease of K as we go to smaller 
smoothing scales similar to that observed in other cross-code 
cross-correlation coefficients. 

One should of course bear in mind that the shape and 
accuracy of the PM force in AP 3 M and ART codes are some- 
what different at small (~ 1 — 4 grid cells) scales. In the 
AP 3 M, the PM force is shaped using modification of the 
Green functions (Couchman 1991) which is controlled by a 
special softening parameter epM- This procedure consider- 
ably increases accuracy of the force (down to < 5% in the 
case of £pm ~ 3.5 used in our PM1 and PM2 simulations) 
in the force, at the expense of making the force somewhat 
"softer" . Thus, for example, in PM1 and PM2 runs, the force 
becomes systematically smaller than the Newtonian value 
at separations <3 grid cells (15% smaller at r = 2 cells 
and 70% smaller at r — 1 cell). The Green functions in the 
PM solver in the ART code are not modified to reduce the 
errors. This results in the force which is Newtonian on aver- 
age down to the scale of one grid cell (the force then falls off 
sharply at smaller separations). At the same time, the errors 
at small separations are considerably higher (see Gelb 1992; 
Kravtsov et al. 1997). At separation of 1 grid cell, the errors 
may reach ~ 50% (albeit in small number of particle config- 
urations). About 10 — 20% of that error is due to the cubical 
shape of the particles assumed in the PM algorithm, while 
the remaining error arises from numerical differentiation of 
the potential. These errors, however, are not systematically 
positive or negative but are scattered more or less evenly 
around zero. This means that particle trajectories can be 
integrated stably down to separations of 1-2 grid cells, as 
was demonstrated in Kravtsov et al. (1997). 

Although it is important to keep the differences in force 
shape and accuracy in mind, we think that they are not 
the main cause of poor cross-correlations. The halo density 
profiles in different PM runs are in good agreement at all 
resolved scales and thus differences in internal density distri- 
butions cannot explain low cross-code K. At the same time, 
visual comparisons of halo positions show small shifts that 
are most naturally explained by the different phase errors 
accumulated in different codes. 

This calls into question the usefulness of cross- 
correlation coefficient for studies of resolution effects (Splin- 



ter et al. 1998), unless the study is done within the same 
numerical code. 

This conclusion can, in fact, be drawn from results of 
Splinter et al. (1998): when the simulation evolves into the 
highly nonlinear stage the cross-correlation within the same 
code is much better than cross-correlation between runs of 
similar resolution but simulated using different codes. For 
example their Table 4 shows that on a 128 3 grid coefficient 
K for the TREE-code runs with the force resolutions of 
e = 0.0625 and e = 0.25 (in the units of mean interparti- 
cle separation) is 0.87, while the cross-correlation coefficient 
between the TREE e = 0.0625 and AP 3 M e = 0.25 runs is 
only 0.67. This is similar to the value of 0.71 which we ob- 
tained for K on the same grid size for ART1 and AP 3 M runs. 
Moreover, Figure 9 in Splinter et al. (1998) also agrees with 
our conclusion. The phase errors demonstrated in this figure 
increase towards smaller scales which explains the decrease 
of cross-correlation coefficient for finer grids. Moreover, the 
figure also shows that cross-correlation within single code 
type is always good regardless of the mass or force reso- 
lution. The largest phase differences are observed between 
runs simulated using different codes. 

The differences in time integration schemes are not 
the only possible sources of small-scale differences in den- 
sity fields. For example, the last two rows of Table [| 
show that cross-correlation between high-resolution ART 
and AP 3 M runs and low-resolution PM run is poor regard- 
less of whether the integration scheme was the same (AP 3 M1 
<-> PM1) or different (ART1 <-> PM1). In this case, the dif- 
ferences in the small-scale details of density fields are due to 
the vastly different force resolutions of the simulations. As 
was noted above, the low resolution of the PM simulation 
(234/i _1 kpc) precludes collapse of any halos with the size 
smaller than ~ 1 — 2 grid cells (see Fig. gj). In the locations of 
small-size halos 8 is very high for ART and AP 3 M runs but 
is much lower (because there are no halos) in the PM run; 
hence the considerably lower cross-correlation coefficient. 

Given the above considerations, our interpretation of 
our results and the results of Splinter et al. (1998) is 
markedly different than the interpretation by the authors 
of the latter study. These authors interpret the differences 
between the low- and high-resolution runs as an erroneous 
evolution in the latter, whereas our interpretation (obvious 
from Fig. jjj) is that these differences are due to the fact that 
high-density small-scale structures such as halos do not col- 
lapse in low-resolution runs. The differences between high- 
resolution runs are interpreted as the phase errors leading 
to small shifts to the locations of small-scale structures, as 
discussed above and as observed in other studies (Frenk et 
al. 1999). 

The origin of these phase errors is the dynamical in- 
stability of particle trajectories in the high-density regions. 
As is well known, the trajectories in the virialized systems 
tend to be chaotic and any small differences existing at any 
time moment will tend to grow very fast with time. The 
divergence can thus be expected to be more important in 
nonlinear regions and this explains the larger phase errors 
at smaller scales. The differences may be caused by the dif- 
ference in the force calculation, errors introduced by numer- 
ical time integration, or simply by different roundoff errors. 
The resulting phase errors are cumulative and thus will grow 
with time. 
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Figure 4. Deviation of particle coordinates at z = in two sim- 
ulations (started with identical initial particle distribution) de- 
pending on the local overdensity (estimated using 128 3 -cell grid) 
at the location of particle in the first simulation. Only 10% of 
all particles are shown. In the upper panel of Fig. ^ we compare 
runs simulated using the same code but different integration steps 
(left) and different force resolution (right). In the lower panel we 
compare runs simulated using different codes. 



Unfortunately, it is impossible to tell which code has 
the smallest phase errors, because we do not know how the 
phases should evolve in the high-density regions. However, 
this is probably not the point. Phase errors of this kind 
will be very difficult to get rid of, because even in the case 
of infinitely good mass, force, and time resolutions, there 
will always be roundoff errors, which will behave differently 
in different codes and thus will tend to grow differences in 
phases. Luckily, almost all popular statistics used in cos- 
mological analyses are not sensitive to phases and therefore 
results are not affected by this problem. Moreover, it is clear 
that at scales > l/i _1 Mpc the errors in phases become neg- 
ligible (different runs cross-correlate perfectly) and there- 
fore even phase-sensitive analyses should not be affected if 
restricted to large scales. However, it is clear that the exis- 
tence of such errors should be kept in mind when analyzing 
or comparing cosmological simulations. 



3.4 Particle trajectories 

As we have discussed in the previous section, small nu- 
merical errors tend to grow and lead to deviations of the 
particle trajectories in nonlinear regions. Nevertheless, it is 
clear that the maximum deviations should be approximately 
equal to the size of a typical halo. Although particle trajec- 
tories can deviate, they are expected to stay bound to the 
parent halo. There is an additional deviation of the order of 
~ 100 — 300ft _1 kpc in the positions of the halos themselves 
(due to the phase errors), but this is also of the order of halo 
size. All in all, we can expect a scatter in the positions of 
the same particles in different simulations not larger than 
the size of the largest systems formed: ~ 1 — 2/i _1 Mpc. 
To compare the particle trajectories in our simula- 
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Figure 5. The same as in Fig. ^ but for the runs AP 3 M1 and 
AP 3 M3. In this case all particles with | Ar > 2h^ 1 Mpc are 
shown (these particles constitute rs 1.4% of the total number of 
particles). The overdensity in the upper panel is estimated at the 
location of the particle in the AP 3 M1 run, while in the bottom 
panel it is estimated for the corresponding particle in the AP 3 M3 
run. Note that counterpart particles in AP 3 M3 tend to be located 
in low-density regions. 



tions, we have calculated the deviations of the coordinates 

Ar = |^l)(fc) _^2)(fc)| ; where ^i)(k) ig the position of the 

k-th particle in the i-th simulation. 

In Fig. [l| we have plotted Ar as a function of the local 
overdensity at the position of the particle for 10% of 
particles randomly selected from the total number of parti- 
cles. In the upper panel of Fig. ^| we compare runs simulated 
using the same code but different integration steps (left) and 
different force resolution (right) . In the lower panels we com- 
pare runs simulated using different codes. 

A quick look at Fig. ^ shows that the scatter in parti- 
cle positions is substantial but in most runs it is contained 
within < 2h^ 1 Mpc, the approximate size of the largest sys- 
tem seen in Fig. |l[ The scatter for most particles in cross- 
code differences shown in the bottom row is somewhat larger 
due to the larger differences in halo positions. As mentioned 
above, this difference adds to the difference in particle posi- 
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Figure 6. Trajectory of two particles in AP 3 M simulation for different force resolutions (At=const). The figure demonstrates the presence 
of two-body scattering in all of our AP 3 M simulations. 
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Figure 7. Acceleration (_in arbitrary units) of particle denoted 
with filled circles in Fig. H. Dotted and dashed curves lowered by 
a constant offset. 



tions within halos. However, the scatter is even larger when 
AP 3 M1 and AP 3 M3 are compared. In this case, about 1.4% 
of particles are separated by more than 2/i -1 Mpc. We did 
not find such outliers when comparing ART runs or 
runs with lower resolution amongst themselves. The compar- 
ison of AP M5 and ART1 shows that the scatter is much 
smaller. 

Figure ^ shows all the particles with separations | Ar > 
IbT 1 Mpc in the AP 3 M1 and AP 3 M3 runs. The overdensity 
in the upper panel is estimated at the location of the particle 
m the AP 3 M1 run, while in the bottom panel it is estimated 
for the corresponding particle in the AP 3 M3 run. It is clear 
that counterpart particles in the AP 3 M3 tend to be located 
in low-density regions, whereas the corresponding particles 
m the AP 3 M1 run are located in a wide range of environ- 
ments. 

The fact that these large deviations occur preferentially 
in the highest resolution runs, immediately raises suspicion 
that they are due to two-body scattering. Indeed, when we 
analyzed the trajectories of the deviant pairs, we found ob- 
vious scattering events. In Fig. ^ we present an example of 
such two-body scattering. In this figure the force resolution 



increases (softening decreases) from right to left. For clarity 
only every second integration step is shown. 

The region shown in Fig. ^ is smaller than the grid 
cell size, so there is no interaction at all in the PM simula- 
tion. The two particles move in the mean potential of the 
other particles. The same happens in the ART simulations. 
On the contrary, due to the high force resolution these two 
particles interact and approach in the case of AP 3 M 1-4 
runs. In the case of the highest force resolution this leads 
to an interaction which, due to the insufficiently small time 
step, leads to a violation of energy conservation. The par- 
ticles approach very closely, feel a large force, undergo in 
this moment a huge acceleration and move away in opposite 
directions with high velocities. At the next integration step 
the particles are already too far to feel substantial two-body 
interaction. Therefore, they move with almost constant high 
velocity in opposite directions. The velocity is about a fac- 
tor of 10 higher than the initial velocity, i.e. the total kinetic 
energy of the system increased during the interaction by a 
factor of 100. Such pairs of particles are located at large Ar 
in scatter plots shown in Figs. ^ & ^| We have run an addi- 
tional simulation with the parameters identical to those of 
the AP 3 M3 run but with a much larger number of time steps 
(48,000 steps in total). The outlying particles dissappear and 
the plots look similar to the AP 3 M5-ART1 plot. This means 
that although the scattering is still present, there is no vi- 
olation of energy conservation in the smaller step run and 
therefore there are no high-velocity streaming particles. 

Fig |B| indicates that scattered particles attain large ve- 
locities and move out of high-density regions. Most of the 
simulation volume is low-density, so it is not surprising that 
we find these streaming particles preferentially in the low- 
density regions. Their counterparts in the AP 3 M1 runs, on 
the other hand, are located in a wide range of environments. 
This implies either that the scattered particles were ejected 
from high-density regions, or that these particles attained 
their excess energy in two-body encounters in low-density 
regions and simply did not participate in the gravitational 
collapse due to their high velocities. 

Fig. shows the acceleration of the particle denoted 
by filled circles in the three AP 3 M plots in Fig. [| This 
figure illustrates the spike in the particle acceleration during 
the scattering event. It also shows that the higher the force 
resolution, the larger the acceleration. 

Such collisions are possible if the force resolution is in- 
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dependent of the local particle density. This problem does 
not arise in the ART code because the resolution is only 
increased in the regions of high local particle density. We 
have checked all particles beyond separation Ar > 2/i _1 Mpc 
within the both ART runs and could not find any event com- 
parable to the collision in the AP 3 M runs. 

The conditions for two-body scattering can be esti- 
mated by noting that strong scattering occurs when the 
potential energy of two-body interaction is equal to ki- 
netic energy of the interacting particles. This gives the scale 
s — 2Gm/v 2 or 



8.61xl0~ 2 /i _1 kpc 



lO 8 /!" 1 M J \ 100 km/s 



.(3) 



For the simulations presented here (rn p — 3.55 x 10 9 h 1 Mq ) 
this scale is 



s = 3.05 x u 10 o h 1 kpc, 



(4) 



where Wioo is velocity in units of 100 km/s. In terms of the 
mean interparticle separation, d = n -1 ^ 3 , this gives 



s/d = 0.013 x v{ 



(•») 



Two-body scattering occurs if the force resolution, e, is 
less than the scale s and if s is much smaller than the local 
interparticle separation di oc . The latter for these simulations 
is d loc = 234.38(1 +<5)~ 1/3 fe~ 1 kpc. 

The above equations show that scattering is possible in 
the AP 3 M runs 1 through 4. This in agreement with results 
presented in Figures ^ and ^. 
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Figure 8. Correlation function of dark matter particles. Note 
that the range of correlation amplitudes is different in the inset 
panel. 



sufficiently small time step, simulate the 2-point correlation 
function incorrectly at scales smaller than the mean inter- 
particle separation. The agreement between high-resolution 
and PM simulations of the same mass resolution always 
agree at the scales resolved in the PM runs. 



3.5 2-point correlation function 

In Fig. [s] we show the correlation function for the dark 
matter distribution down to the scale of 5/i _1 kpc, which 
is close to the force resolution of all our high-resolution 
simulations. The correlation function in runs AP 3 M1 and 
ART2 are similar to those of AP 3 M5 and ART1 respec- 
tively and are not shown for clarity. We can see that the 
AP 3 M5 and the ART1 runs agree to < 10% over the whole 
range of scales. The correlation amplitudes of runs AP 3 M 2- 
4, however, are systematically lower at r<50 — GOh^ 1 kpc 
(i.e., the scale corresponding to ~ 15 — 20 resolutions), with 
the AP 3 M3 run exhibiting the lowest amplitude. At scales 

< 30/i _1 kpc the deviations from the ART1 and the AP 3 M5 
runs are « 100 — 200%. We attribute these deviations to the 
numerical effects discussed in § 5. The fact that the AP 3 M2 
correlation amplitude deviates less than that of the AP 3 M3 
run, indicates that the effect is very sensitive to the force 
resolution. 

The correlation function of the PM simulation deviates 
strongly on small scales. However, the bend coincides with 
the force resolution (~ 234/i~ 1 kpc) of this run. At scales 
smaller than resolution, we can expect an incorrect corre- 
lation amplitude because waves of wavelength smaller than 
the resolution do not grow at the correct rate in these runs. 

This result agrees with the correlation function com- 
parison done by Cohn et al. (1999), where agreement of 

< 10% was found between the correlation functions from 
the larger 256 3 -particle ART, AP 3 M, and PM simulations 
at all resolved scales. There is no evidence therefore, that 
high-resolution simulations, given that they are done with 



4 PROPERTIES OF HALOS 
4.1 Identification of Halos 

To compare properties of the halos and their distribution 
in different simulations, we identify DM halos using the the 
friends-of-friends (FOF) algorithm. The algorithm identifies 
clumps in which all particles have neighbors with distances 
smaller than times the mean inter particle separation, 
ru < 11 x I. This halo finding algorithm does not assume 
a special geometry for the identified objects. A drawback is 
that it only uses particle positions and therefore can identidy 
spurious unbound clumps at low halo masses. 

The mean overdensity of a particle cluster is related to 
the linking length used to identify it. In an Einstein-de-Sitter 
universe (simulated in our runs) virialized halos have an 
overdensity of 5 V [ T m 178, which corresponds approximately 
to the linking length of 11 — 0.2. A reduction of the linking 
length by a factor of 2 roughly corresponds to an increase of 
the overdensity by a factor of 8. Using smaller linking lengths 
we can study the substructure of the DM halos. Table |§] lists 
the number of halos identified using different linking lengths 
in different runs for the two lower limits on the number of 
particles in a halo: 25 (columns 2 — 4) and 50 (columns 5 — 7). 

The table shows that there are differences of ~ 15 — 50% 
between high-resolution runs. The differences are present 
not only at low linking lengths, but even at the "virial" link- 
ing length, // = 0.2. These differences are partly due to the 
nature of the FOF algorithm: small differences in the par- 
ticle configurations may lead to an identification of a single 
halo in one simulation and to the identification of two or 
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Table 3. Number of DM halos identified using the FOF algorithm 
as a function of the linking length (in units of mean interparticle 
separation). Only halos containing more than 25 (columns 2 — 4) 
and 50 (columns 5 — 7) particles have been counted. 
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linking length: 


0.2 


0.1 


0.05 


0.2 


0.1 


0.05 


AP 3 M5 


388 


292 


162 


206 


152 


85 


AP 3 M1 


377 


264 


160 


204 


149 


84 


AP 3 M4 


340 


221 


110 


187 


120 


65 


AP 3 M2 


365 


251 


122 


188 


132 


74 


AP 3 M3 


332 


215 


87 


173 


117 


55 


PM1 


214 


150 


50 


134 


74 


19 


PM2 


208 


155 


49 


129 


66 


18 


ART1 


359 


291 


161 


199 


148 


89 


ART2 


383 


297 


169 


206 


156 


87 



more halos in another simulation. Nevertheless, some sys- 
tematic differences are also apparent. PM simulations have 
almost half as many halos due to the absence of small-mass 
halos (see Fig. [l]) that do not collapse (or do not survive) 
due to the poor force resolution of these runs. 

Also, while AP 3 M1, AP 3 M5, ART1, and ART2 runs 
agree reasonably among themselves, the number of halos in 
runs AP 3 M2, AP 3 M3, and AP 3 M4 is systematically lower. 
In this case the differences seem to be counter-intuitive: 
the number of halos is smaller in higher force resolution 
runs. These differences persist at all linking lengths indi- 
cating that there are differences in substructure as well as 
in the number of isolated halos. Moreover, the differences 
persist even for a larger cut in the number of particles. It 
has been noted in previous studies (e.g., van Kampen 1995; 
Moore, Katz & Lake 1996) that particle evaporation due to 
two-body scattering (see, for example, Binney & Tremaine 
1987) can be important for halos of < 30 particles. For such 
halos the evaporation time-scale, especially in the presence 
of strong tidal fields, can be comparable to or less than the 
Hubble time. However, for halos containing > 50 particles, 
evaporation should be negligible. Nevertheless, the trend 
with resolution seen in Table |^ does suggest that two-body 
evaporation is the process responsible for the differences. 

The most likely explanation of these results is, in our 
opinion, the accuracy of the time integration. The estimates 
of the evaporation time-scale are done assuming no errors 
in energy exchange between particles in a scattering event. 
Our analysis of such events in our simulations, on the other 
hand, shows that with the time step of the AP 3 M 2-4 runs, 
there is severe energy conservation violation during scatter- 
ings. The particles attain much larger velocities than they 
should have if the integration were perfect. For example, in 
the scattering event shown in Fig. m the final kinetic energy 
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Figure 9. Mass function for some simulations and two linking 
lengths. 

of the two particles is 100 times larger than their initial 
kinetic energy. Therefore, a time step that is insufficiently 
small to properly handle two-body scattering may exacer- 
bate the process of evaporation and result in much shorter 
than predicted evaporation time-scales, even for halos con- 
taining relatively large numbers of particles. The differences 
between AP 3 M1 run and runs AP 3 M2 and 4 indicate that 
this effect is very sensitive to both the force resolution and 
to the time step of the simulation. 

Another possible explanation is that halos do not evap- 
orate but the particles are heated due to non-conservation 
of energy which makes halos "puffier" . Such halos would be 
more prone to destruction by tides in high-density regions. 

4.2 Mass Function 

In Fig. ^ we show the mass functions of all halos identified 
with various linking lengths at z — 0. 

As can be seen in the left panel of Fig. [], the halo mass 
function in the PM run is biased towards large masses: the 
number of objects of mass >2 x 10 12 /i _1 M in the PM 
run agrees well with the corresponding number in the high- 
resolution runs, while at lower masses the number of ha- 
los is strongly underestimated. Mass 2 x 10 12 /i _1 Mq cor- 
responds to the virial radius (at z = 0, <5 v ir = 178) of 
R v ir ~ 213/i _1 kpc, i.e. very close to the force resolution 
of the PM runs (234/i _1 kpc). The conclusion is that the 
PM runs fail to reproduce the correct abundances of halos 
with the virial radius less than about two grid cell sizes. 

At smaller linking length, 11 = 0.05, (right panel of 
Fig. ^ the PM run severely underestimates the halo abun- 
dances. We attribute this also to the poor force resolution 
of the run. The poor resolution prevents the formation of 
dense cores in the inner regions of collapsed halos and ha- 
los that do not reach overdensity of » 11,000 (overdensity 
of objects identified with 11 = 0.05) will be missed. Even if 
the central density of some halos reaches this value, the ha- 
los are still considerably less dense than their counterparts 
in the high-resolution runs and are therefore more suscepti- 
ble to destruction by the tidal fields in high-density regions 
(Moore et al. 1996; Klypin et al. 1999). 

The mass functions of the high-resolution runs agree to 
~ 30% for isolated halos of overdensity 5 = 178 (11 — 0.2). 
The AP 3 M runs have a larger number of identified halos at 
masses < 5x 10 10 /i _1 Mq (i.e., halos containg < 15 particles) 
than the ART runs. This difference is due to smaller number 
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Figure 11. Correlation function of halos for some simulations 
and two linking lengths. A mass cut of 25 particles was used to 
construct halo catalogs. 



of small halos in low-density regions in the ART runs dis- 
cussed in § 3.2. The mass functions of the AP 3 M1, AP 3 M5, 
and the ART runs agree well at masses >10 n /i _1 Mq for 
both // = 0.2 and 11 = 0.05 (overdensity of 178 and 13000, re- 
spectively), indicating that both runs produced similar pop- 
ulations of halos with similar central densities. The mass 
functions of the AP 3 M 2-4 runs are similar for ll = 0.2, but 
show differences for ll — 0.05. Thus, for example, the abun- 
dance of halos of mass ~ 10 11 - 10 12 /i -1 M (~ 30 - 300 
particles) in the AP 3 M3 run is underestimated by a factor 
of w 1.5 - 2. The mass functions of the AP 3 M2 and AP 3 M4 
runs lie in between those of the AP 3 M3 and AP 3 M5 runs. 

The fact that differences are present at small linking 
length indicates differences in the high-density regions. This 
may be due to the generally lower inner densities of halos 
and/or to the destruction of "heated" sattelites discussed 
above. 



4.3 Halo correlation function 

Figure [ll] shows the 2-point correlation function of identi- 
fied DM halos. There is good agreement between the cor- 
relation functions of isolated virialized (ll = 0.2) halos in 
high-resolution runs. Similar to the dark matter correlation 



function, the agreement is better than 10%. The agreement 
between the AP 3 M1, AP 3 M5, and the two ART runs does 
not break even at higher overdensities (ll = 0.05), which 
indicates that these runs produced similar small-scale sub- 
structures in the high-density regions within isolated halos. 
We do not find any differences of the type seen in the DM 
correlation function (§ 3.5) between AP 3 M runs at = 0.2. 
For halos identified using ll — 0.05 some differences are ob- 
served, but these are comparable to the poisson errors and 
are therefore inconclusive. 

On the other hand, halos in the PM simulation exhibit 
higher correlations than halos in the high-resolution runs. 
As discussed in the previous section, the halo mass function 
in the PM run is biased toward high masses. The higher 
amplitude of the correlation function can thus be explained 
by the mass-dependent bias: higher mass halos are clustered 
more strongly. 



4.4 Velocity Dispersion vs. Mass 

Figure |l(j shows the correlation of velocity dispersion and 
mass for one of the ART, AP 3 M, and PM simulations. A 
correlation ov oc M 1//3 is expected for virialized halos. For 
the ART1 and AP 3 M5 simulations, the best fit slope of the 
correlation for the 50 most massive halos is 0.33 (the corre- 
lation for the rest of the AP 3 M runs is similar). The velocity 
dispersion of low-mass halos scatters around this value. For 
comparison, a line of slope 1/3 is included in all three panels. 
It should be mentioned here that at the low mass end the 
FOF groups contains only a few particles so that o v is not 
well determined because the error due to unbound particles 
accidentally linked by the FOF algorithm may be very high. 

The halos in the PM run deviate from the expected 
slope at masses < 10 13 /i -1 Mq. The virial radii of the halos 
of these masses are <364/i _1 kpc. Their size is thus < 3 force 
resolutions across. Therefore, the potential and the internal 
dynamics of the particles in these halos are underestimated 
by the PM code leading to a steeper slope of the ov — M 
relation. 
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Figure 12. Density profiles of halos in different runs. 



4.5 Density Profiles 

In this and the next two sections we present halo-to-halo 
comparison of individual halos in different simulations. In 
this section we will compare the density profiles of DM ha- 
los. The density distribution of hierarchically formed halos 
is currently a subject of debate (see, e.g., Navarro, Frenk & 
White 1997; Kravtsov et al. 1998; Moore et al. 1999) and 
study of the resolution effects and cross-code comparisons 
are therefore very important. In figure [l2| we present the 
density profiles of four of the most massive halos in our sim- 
ulations. We have not shown the profile of the most massive 
halo because it appears to have undergone a recent major 



merger and is not very relaxed. In this figure, we present only 
profiles of halos in the high-resolution runs. Not surprisingly, 
the inner density of the PM halos is much smaller than in 
the high-resolution runs and their profiles deviate strongly 
from the profiles of high-resolution halos at the scales shown 
in Fig. n3. We do not show the PM profiles for clarity. 

A glance at Fig. ^ shows that all profiles agree well at 
r > SOh^ 1 kpc. This scales is about eight times smaller than 
the mean interparticle separation. Thus, despite the very 
different resolutions, time steps, and numerical techniques 
used for the simulations, the convergence is observed at a 
scale much lower than the mean interparticle separation, ar- 
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Figure 13. Cumulative number of particles within radius r from the centers of four halos shown in the previous figure. 



gued by Splinter et al. (1998) to be the smallest trustworthy 
scale. At smaller scales the profiles become more noisy due 
to poorer particle statistics (see Fig. |l3| ). 

Nevertheless, there are systematic differences between 
the runs. The profiles in two ART runs are identical within 
the errors indicating convergence (we have run an additional 
simulation with time steps twice smaller than those in the 
ART1 finding no difference in the density profiles). Among 
the AP 3 M runs, the profiles of the AP 3 M1 and AP 3 M5 are 
closer to the density profiles of the ART halos than the 
rest. The AP 3 M2, AP 3 M3, and AP 3 M4, despite the higher 
force resolution, exhibit lower densities in the halo cores, the 



AP 3 M3 and AP 3 M4 runs being the most deviant. These dif- 
ferences can be seen more clearly in Fig. [l3|, where we plot 
the cumulative number of particles (i.e., mass) within radius 
r for the halos shown in Fig. The differences between 
AP 3 M3 and AP 3 M4 and the rest of the runs are apparent 
up to the radii containing ~ 1000 particles. 

These results can be interpreted, if we examine the 
trend of the central density as a function of the ratio of 
the number of time steps to the dynamic range of the simu- 
lations (see Table |l]) shown in Table [|. The ratio is smaller 
when either the number of steps is smaller or the force res- 
olution is higher. Table H shows that agreement in density 
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Table 4. The ratio of number of time steps to the dynamic range 
of simulations for the high-resolution runs. The number of time 
steps for the ART runs correspond to the time step on the highest 
refinement level. 



Run 


N step /dyn. range 


AP 3 M4 


0.469 


AP 3 M3 


0.702 


AP 3 M2 


0.938 


AP 3 M1 


1.875 


AP 3 M5 


3.750 


ART2 


2.578 


ART2 


5.156 



profiles is observed when this ratio is > 2. This suggests that 
for a fixed number of time steps, there should be a limit on 
the force resolution. Conversely, for a given force resolution, 
there is a lower limit on the required number of time steps. 
The exact requirements would probably depend on the code 
type and the integration scheme. For the AP 3 M code our 
results suggest that the ratio of the number of time steps to 
the dynamic range should be no less than one. It is interest- 
ing that the deviations in the density profiles are similar to 
and are observed at the same scales as the deviations in the 
DM correlation function (Fig. []) suggesting that the corre- 
lation function is sensitive to the central density distribution 
of dark matter halos. 

These results are indicative of the sensitivity of the cen- 
tral density distribution in halos to parameters of the numer- 
ical simulation. However, due to the limited mass resolution 
of our test runs, they do not shed light on the density profile 
debates. The profiles of the ART halos agree well with those 
of the AP 3 M halos, if the latter are simulated with a suffi- 
ciently large number of time steps. But debated differences 
are at scales of 2 0.01.Rvir, which are not resolved in these 
simulations. We are currently carrying out a more detailed, 
higher-resolution study to clarify the issue. 

4.6 Shared Particles 

We have argued above that particle trajectories may diverge 
in high-density regions due to their instability to small inte- 
gration errors. However, despite the instability of its trajec- 
tory, we can expect that a bound particle will stay bound 
to its parent halo. To this end we compare the particle con- 
tent of individual halos in different runs. We have identified 
the counterpart halos in different runs as systems that have 
the largest number of common (or shared) particles. This 
method is superior to coordinate based methods because 
due to the phase errors (see § 3.3) one does not expect to 
find the halos at the exact same position in different runs 
(this is especially not for small groups of particles). 

Fig. [14] shows the percentage of mass that can be found 
in particle groups that coincide in a specific amount of 
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Figure 15. Correlation of the velocity dispersions of the same 
halos identified in two different simulations. Only halos with more 
than 50 particles are taken into account 

shared particles. The figure shows that in all runs most ha- 
los have more than 80% of their particles in common. The 
distributions peak at around 85-90% for all of the compar- 
isons, except ART1 vs. ART2 which peaks at « 90 - 95%. 
Even though a direct comparison of particle positions within 
halos is not a very useful way of comparing different runs, 
this result shows that the particle content of halos is rather 
similar. 



4.7 Correlation of Velocity Dispersions 

In Fig. |l5| we compare the velocity dispersions of halos iden- 
tified in different simulations by the FOF algorithm assum- 
ing a minimum particle number of 50 per halo. The figure 
shows that the velocity dispersions agree reasonably well 
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with the overall scatter of ~ 50 km s . The very few out- 
liers are halos of very different masses identified as the same 
halo by the FOF algorithm. Small differences in particle dis- 
tribution may result in the identification of a binary system 
in one simulation and only a single halo in another simula- 
tion, leading to the large differences in velocity dispersions. 

The differences between velocity dispersions appear to 
be independent of halo mass, and although they reach 
~ 50% for the lowest mass halos, there are no obvious sys- 
tematic differences between different simulations. 



5 DISCUSSION AND CONCLUSIONS 

We have presented results of a study of resolution effects in 
dissipationless simulations. As we noted in the introduction, 
an additional goal of this study was to compare simulations 
done using two different high-resolution TV-body codes: the 
AP 3 M and the ART. Our results indicate that both codes 
produce very similar results at all scales resolved in the pre- 
sented simulations, given the force resolution and time step 
are such that convergence is reached within the code type. 
Our results indicate that numerical effects may be compli- 
cated due to a combined effect of mass and force resolution 
and inaccuracies of the time integration scheme. The precise 
magnitude of the effects depends on the numerical parame- 
ters used and the considered statistics or property of particle 
distribution. 

Particles in dissipationless cosmological simulations are 
supposed to represent elements of the dark matter distribu- 
tion in phase space The smaller the particle number, the 
larger the volume associated with each of the particles. Dur- 
ing the course of evolution, according to Liouville's theorem, 
the phase-space volume of each element should be preserved. 
Its shape, however, will be changing. Correspondingly, the 
eulerian space volume of these elements may shrink (for par- 
ticles in an overdense collapsing region of space) or expand 
(for particles in underdense regions). Regardless of its initial 
shape, each element can be stretched due to the anisotropic 
nature of the gravitational collapse in cosmological models. 

Usually, none of these effects is modelled in cosmologi- 
cal TV-body simulations. The gravitational field of each par- 
ticle is assumed to be roughly isotropic and therefore the 
effects of the volume stretch cannot be modelled. These ef- 
fects are not addressed in the present study, because they 
can only be studied by widely varying the number of par- 
ticles at a fixed force resolution and for a fixed time step^J 
while we have varied the force resolution and time step keep- 
ing the number of particles fixed. A large number of parti- 
cles corresponds to the smaller phase-space associated with 
each particle and the symmetric particle approximation is 
more accurate. Convergence studies of the halo density pro- 
files indicate that these effects are small (at least at radii 
> 0.02 — 0.05-Rvir)- However, this does not mean that they 
cannot be important for other statistics. 



t Unless particles are supposed to represent individual galaxies, 
which is virtually never the case in modern simulations, 
t Due to the phase errors discussed below, such a study should 
also be carried out using the same numerical code. 



Nevertheless, we can study other kinds of resolution ef- 
fects. Normally, softening of interparticle gravitational force 
should be approximately equal to the spatial size of the 
phase-space volume element associated with each particle. 
If the softening is much smaller than this size, the volume 
elements will behave like particles and two-body scattering 
is possible. This was indeed observed in some of our high- 
resolution simulations. This scattering contradicts the col- 
lisionless nature of the modelled dark matter and is thus 
undesirable. 

For the mass resolution of our simulations (~ 3.55 x 
10 9 /i _1 Mq), scattering was observed in simulations with 
uniform force resolution of <3/i _1 kpc and disappears for 
larger values of resolution. While we find strong scattering 
for ~ 1.4% of the particles, many more may have suffered 
weaker scattering events during the course of simulation. 
Indeed, we find indirect evidence that scattering has sub- 
stantial effect on the particle distribution. In particular, we 
find that it may noticeably affects the small-scale ampli- 
tude of the 2-point correlation function of dark matter, the 
abundance and mass function of dark matter halos, and the 
central density distribution of halos. The effect is amplified 
strongly for simulations in which larger time steps are used. 
This is because for larger time steps the scattering events 
severely violate energy conservation. The magnitude of the 
violation is very sensitive to the time step: the amount of 
momentum gained by a particle during a scattering event 
is oc gAt, where g is the acceleration. Therefore, for the 
same force resolution (which means that the same g can be 
achieved), the momentum gain will be proportional to At. 

The 2-point correlation function of matter and the cen- 
tral density distribution in halos in our runs are affected by 
these effects at scales as large as ~ 15 — 20 force soften- 
ings. Also, the abundance of satellite halos in high-density 
regions appears to be affected as well. We think that these 
effects may be due to the following two phenomena. First, 
if the time step of the simulation is too large for a given 
force resolution, particle trajectories are not integrated ac- 
curately in the highest density regions, where the gradient 
of potential is the highest. In some sense, particles scatter 
off the central density peak, and may gain energy when in- 
tegrated with large time step (see arguments above). The 
number of simulation time steps per deflection on a scale R 
can be estimates as 

N stcp = ~~\7 « 7 - 51 x W~ 4 RK?e P vToo, (6) 

where R is in kpc, At and TV s °e P is the time step and the 
total number of time steps of the simulations, and wioo is 
particle velocity in units of 100 km s -1 . We have assumed 
the Einstein-de Sitter universe with the Hubble constant of 
Ho — 50 km s _1 Mpc -1 . For high- velocity particles stream- 
ing through the centers of massive halos, there may be just 
a few time steps to integrate the part of trajectory where 
strong changes in acceleration occur. As illustrated in Fig. [| 
this may not be sufficient to ensure energy preservation and 
may result in an energy gain by particles. This leads to arti- 
ficial heating and lowers the central density because having 
acquired energy, particles are not as likely to enter the cen- 
tral region of halo. 

The second phenomenon is due to the "graininess" of 
the potential. Particles in the high-density regions may feel 
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the discreteness of the density field and suffer scattering. We 
do find evidence for scattering in the high-density regions in 
our simulations (see § 3.4). Here, again, the effect may be 
amplified strongly by the incorrect integration of such scat- 
tering. Indeed, run AP 3 M1 (N^ p = 8000) performs much 
better than the run AP 3 M4 (JV* t e p = 2000), although both 
runs have the same mass and force resolution. 

The two effects described above may operate in com- 
bination, although the first effect does not depend on the 
mass resolution. The second effect should be eliminated for 
higher mass resolution. Both lead to the artificial heating of 
particles thereby lowering the central density of halos and 
possibly ejecting the particles altogether in some cases. Vi- 
sual comparisons of halos in the AP 3 M3 (the run which per- 
formed the worst) and the AP 3 M5 runs shows that AP 3 M3 
halos appear "puffier" and more extended than the same 
halos in the AP 3 M5 or ART runs. Puffier halos may be de- 
stroyed more easily by tides in high-density regions which 
may explain some of the differences seen between the mass 
functions of halos in different runs. 

Our results show that in constant time step high- 
resolution simulations the total number of time steps must 
be rather high to ensure good energy conservation. This re- 
quirement can become computationally prohibitive in sim- 
ulations that follow large numbers of particles. In case of 
the AP 3 M code, it would probably be preferrable to use its 
version in the publicly available code "Hydra" (Couchman 
et al. 1995) that uses adaptively varied time step. 

The conditions for scattering, discussed in § 3.4, occur 
if the force softening is smaller than the scale s, which in 
units of mean interparticle separation is 

i «-'« , (is^)"'( i j^) W w 

and is considerably smaller than the local interparticle sepa- 
ration: di oc = (1 + 5) -1 ^ 3 (in units of the mean interparticle 
separation, S is the local particle overdensity) . For our sim- 
ulations s ~ 3h^,q/i _1 kpc, where vioo is particle velocity 
in units of 100 km s _1 . This means that condition s -C di oc 
is satisfied everywhere but in the highest density regions: 
5>W 4 . The conditions for strong scattering occur for the 
slow moving (<100kms _1 ) particles in the AP 3 M runs 
1 — 4. Such slow moving particles are likely to be present 
in the low-density regions and in small-mass halos of ve- 
locity dispersion <r v < 100 — 200 km s -1 . This may explain 
our result that halos of mass < (0.5 - 1) x 10 12 h _1 M 
(<r v ~ 100 — 150 km s _1 , see Fig. lid) appear to be affected 
by scattering. 

One may question the relevance of these results, given 
the small size of the simulations and extremely high force 
resolution. Note, however, that our results would be appli- 
cable (save for the presence of very massive clusters) to any 
256 3 -particle simulation, in which force resolution of consid- 
erably smaller than the scale s is adopted. The simulations 
with the particle mass and dynamic range not very far from 
ours, have already been done. For example, all 256 3 -particle 
simulations of 239. 5/i -1 Mpc box presented by Jenkins et al. 
(1998) satisfy the condition of e < s (where s is estimated 
using the above equation for wioo S: 1-5) and have been run 
using < 1600 time steps. The parameters of the recent "Hub- 
ble volume" simulation (Colberg et al. 1998) also satisfy con- 



ditions for strong scattering: s 1.9/i _1 Mpc, while force 
resolution of the simulations is e = 100/i -1 kpc. The parti- 
cle mass of these simulations is ~ 2 x 10 12 /i -1 Mq, which 
means that individual particles represent galaxies rather 
than a phase-space element. Galaxies form a collisional sys- 
tem so the presence of scattering may be considered as a 
correct model of the evolution of galaxy clustering. More- 
over, the mean interparticle separation in these simulations 
is ~ 2 — 3/i _1 Mpc and thus conditions for scattering may 
only occur in the underdense regions. 

A high-resolution 256 3 -particle run was also done re- 
cently using the ART code (see, for example, Coh'n et 
al. 1999). However, for this simulation (m p = 1.1 x 
10 9 fe _1 Mq; f2o = 0.3) the scale of strong scattering is 
s ~ 0.94hJ~ q/i _1 kpc, while peak resolution is « 4/i -1 kpc. 
The time steps for the particles at the refinement level of 
this resolution corresponds to effective number of time steps 
of ~ 41,000. Therefore, for this simulation the strong scat- 
tering condition is not satisfied. Moreover, a refinement level 
L is introduced in these simulations only if the local over- 
density is higher than 8 = 5 x 2 3 ^ L+1 \ or for the highest 
resolution level L = 6: S » 10 7 . For these overdensities, the 
local interparticle separation is ~ lh, -1 kpc, and two-body 
interactions are thus unlikely. 

To summarize, scattering can be precluded if the choice 
of force resolution is guided by the scale s, which, in turn, 
depends on the mass resolution (particle mass). This con- 
clusion may seem similar to that of Splinter et al. (1998), 
who concluded that force resolution should not be smaller 
than the mean interparticle separation. It is, however, quite 
different in practice: eq. (7) shows that for our box size 
s ~ 1 in units of mean interparticle separation only for 
m p Ri3x W^h' 1 M e . For our mass resolution, force soft- 
ening as small as 5 — 10h~ kpc is justified. This is ~ 25 — 50 
times below the mean interparticle separation. 

We think that the conclusion of Splinter et al. is (at least 
in part) due to the interpretation of poor cross-correlation 
between different simulations on small scales as erroneous 
evolution in high-resolution runs. Our analysis, presented 
in § 3, shows that poor cross-correlation is due to phase 
errors whose major source is cumulative errors due to the 
inaccuracies of time integration. The trajectories of particles 
become chaotic in high-density regions and small differences 
in time integration errors tend to grow quickly. 

For this reason it may prove to be very difficult to get 
rid of this effect by improving the time integration. There- 
fore, one should keep these errors in mind if a phase sen- 
sitive statistic is analyzed. Luckily, most of the commonly 
used statistics are phase-insensitive and are not affected by 
such errors. Moreover, the errors are confined to the small- 
scale high-density regions, and no significant phase errors are 
present in our simulations if the density field is smoothed on 
a scale -J l/i -1 Mpc. 

While this is clearly still an error, it has nothing to do 
with the mass or force resolution and would be present even 
if both were perfect. This point is clearly demonstrated by 
the fact that simulations run using two different implemen- 
tations of the PM code correlate perfectly within the code 
type but cross-correlate rather poorly when cross-code com- 
parisons are made (see § 3.3). Note that in all of these PM 
runs, the force resolution is approximately equal to the mean 
interparticle separation. 
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The main conclusion of our study is that care must be 
taken in the choice of force resolution for simulations. If a 
code with spatially uniform force resolution is used, condi- 
tions for strong two-body scattering may exist if the force 
resolution is smaller than the scale s discussed above. The 
presence of scattering itself may not be important (albeit 
undesirable); the relaxation time for systems, for example, 
may be much longer than the Hubble time (e.g., Hernquist 
& Barnes 1990; Huang, Dubinski & Carlberg 1993). Its ef- 
fects, however, may be greatly amplified if the time step of 
the simulation is not sufficiently small. In this case, severe 
violation of energy conservation occurs during each scatter- 
ing which may lead to artificial injection of energy into the 
system. 
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